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ABSTRACT 


This  paper  is  concerned  with  a  numerical  study  of  the  equation 


Pu  =  3n 


3x3t 


u  (t) 
Kt-s)(-§— 
ux(s) 


u  (s) 
x  )ds 


U2(t) 

X 


where  u(x,t)  is  real  valued  for  x  e  [-1,1]  and  t  e  R  with  the  boundary 
condition 


u  (t) 
a(t-s) (-^ - 

Ux(8) 


- - )dt  =  f(t) 

ux(t) 


at  x  =  ±1 .  This  problem  is  a  model  equation  for  elongation  of  a  thin 
filament  of  a  polymeric  liquid  when  the  force  f  is  applied  at  both  ends. 

The  initial  condition  is  u(x,-")  =  V>(x) .  The  unknown  variable  u(x,t) 
denotes  the  position  of  a  fluid  particle  (in  a  deformed  state  at  time  t) , 
which  is  at  position  x  in  space  in  a  certain  reference  configuration.  In 
this  reference  state  the  filament  is  assumed  to  be  cylindrical.  a(t)  is  a 
memory  kernel,  P  denotes  the  density  of  the  fluid  and  n  the  Newtonian 
contribution  to  the  viscosity.  We  set  up  a  difference  scheme  for  this  problem 
and  show  the  convergence  under  certain  assumptions  on  f  and  we  report 
computations  . 
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SIGNIFICANCE  AND  EXPLANATION 


This  paper  is  concerned  with  the  study  of  a  model  equation  for  the 
elongation  of  a  filament  of  a  polymeric  liquid  when  a  force  is  applied  to  the 
ends  of  the  filament.  Mathematically  the  equation  has  the  form  of  a  nonlinear 
partial  integrodif f erential  equation  for  the  position  of  a  fluid  particle  in 
time  and  space.  However  this  equation  can  be  transformed  to  a  quasilinear 
parabolic  system.  Parameters  of  the  equation  are  the  density  of  the  fluid, 
the  Newtonian  contribution  to  the  viscosity  and  certain  relaxation  constants 
which  represent  the  'memory'  of  the  fluid. 

Our  focus  is  the  numerical  investigation  and  therefore,  a  difference 
method  is  set  up  and  its  convergence  is  proved  under  the  assumption  that  the 
force  which  is  applied  to  the  ends  of  the  filament  is  small  in  an  appropriate 
norm  (this  is  the  case  for  which  an  existence  and  uniqueness  theorem  was 
proved)  .  Computations  are  reported. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC,  and  not  with  the  authors  of  this  report. 


THE  NUMERICAL  SOLUTION  OF  A  QUASILINEAR  PARABOLIC 
EQUATION  ARISING  IN  POLYMER  RHEOLOGY 


P.  Markowich* '  ^  and  M.  Renardy  ' ^ 


1 .  INTRODUCTION 


We  present  a  numerical  study  of  a  model  equation  describing  the 
stretching  of  a  filament  of  a  polymeric  liquid,  when  a  force  f  is  applied  to 
both  ends: 


The  model  is  based  on  the  following  assumptions:  The  polymer  is 
incompressible  and  fulfills  the  "rubber-like  liquid"  constitutive  relation 
[3] .  Moreover,  the  filament  is  thin  and  hence  the  originally  three- 
dimensional  problem  can  be  approximated  (formally)  by  an  equation  which  is 
one- dimensional  in  space. 

Under  these  assumptions,  the  following  model  equation  was  derived  in  [7] . 


(1.1) 
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,  „  «  t  u  (t)  u  (a) 

(1.2)  3i»|^(-l-)+  /  a(t-s)  (— - * - )da  -  f(t) 

Ux  --  u(s)  u2(t) 

at  x  -  11,  t  e  R  X  X 

(1.3)  u(x,t  *  -*)  -  *>(x) 

for  x  e  [-1,1] 

The  subscript  x  denotes  partial  differentiation  with  respect  to  x,  and 

"dot"  denotes  partial  differentiation  with  respect  to  t.  The  unknown 

variable  u(x,t)  is  the  position  of  a  fluid  particle  in  a  "deformed"  state 

(at  time  t),  which  is  at  position  x  in  an  "undeformed"  reference  state.  In 

this  reference  state  the  filament  is  assumed  to  be  cylindrical.  Thus  the 

thickness  of  the  filament  in  the  initial  state  u  *  ^(x)  is  proportional  to 
1 

v>*(x)  * 

We  have  normalized  the  length  scale  such  that  the  filament  has  length 
2  in  the  reference  state.  The  resulting  scaling  factor  is  absorbed  into  Pi 
which  denotes  the  density  of  the  filament  multiplied  by  the  square  of  half  the 
length  in  the  reference  state .  n  denotes  a  Newtonian  contribution  to  the 
viscosity,  and  f  the  force  acting  on  the  ends  divided  by  the  cross-sectional 
area  in  the  reference  state . 

The  memory  kernel  a  has  the  form 

M  -X  t 

(1.4)  a(t)  “  I  K  e  ,  K  ,  X  >  0  . 

i-1  1 

(more  general  kernels  were  considered  in  [7] ) . 

A  problem  closely  related  to  (1.2)  was  investigated  analytically  by 
Lodge,  Me  Leod  and  Nohel  [4]  and  numerically  by  Nevanlinna  [6]  .  They  regard 
(1.2)  as  a  history  value  problems  ux  is  given  for  t  >  0,  and  f(t)  »  0 
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for  t  >  0.  Their  results  hold  for  a  class  of  kernals  a  and  functions 

F(ux(t),  u^fs))  under  the  integral  that  include  those  in  (1.2). 

In  our  current  setting,  (1.2)  was  investigated  both  analytically  and 

numerically  in  [5]  as  a  model  for  the  elongation  of  the  filament,  when 

inertial  forces  are  neglected.  It  was  shown  there  that  if  f  is  small  in  the 

00 

sense  of  an  exponentially  weighted  L  -norm  or  if  f  converges  to  0 
exponentially  as  t  ♦  -*•  and  f(t)  =  0  for  t  >  tQ ,  t  finite,  then  problem 
(1.2)  has  a  unique  solution  globally  in  time  which  satisfies  the  initial 
condition  at  t  =  -°° .  This  solution  lies  again  in  an  exponentially  weighted 

oo 

L  -space.  In  [7],  Renardy  showed  that,  for  the  case  of  small  f,  a  similar 
existence  result  holds  for  the  full  problem  (1.1).  This  equation  was 
transformed  to  a  quasilinear  parabolic  system  that,  together  with  (1.2), 

2 

constitutes  a  Neumann  problem.  For  technical  reasons,  weighted  Sobolev  (L  - 

00 

based)  spaces  were  used  instead  of  weighted  L  -spaces  here. 

For  the  boundary  problem  (1.2),  it  was  also  shown  that  solutions  depend 
continuously  on  n  even  as  n  +  0.  No  result  like  this  is  yet  known  for  the 
full  spatial  problem. 

Although  this  does  not  affect  the  mathematics,  we  wish  to  remark  that 
only  f  >  0  makes  sense  physically,  because  an  attempt  to  compress  the 
filament  would  result  in  buckling. 

In  this  paper,  we  devise  a  finite  difference  method  for  the  numerical 
solution  of  (1.1)  -  (1.3),  and  we  show  the  convergence  of  this  method  for 
small  forces  f  and  the  initial  condition  v>(x)  =  x.  Our  proof  uses  discrete 
analogues  of  the  weighted  Sobolev  spaces  employed  in  the  analytical  theory. 

The  paper  is  organized  as  follows.  Chapter  2  recalls  the  existence 
result  given  in  [7] ,  chapter  3  contains  the  convergence  proof  for  the 
difference  scheme,  and  computations  are  reported  in  chapter  4. 
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2.  ANALYTICAL  THEORY 

We  study  the  problem  (1.1)  -  (1.3)  with  a  kernel  of  the  form  (1.4).  The 
existence  result  we  present  Is  "local"  In  the  sense  that  f  Is  "small"  in  a 
certain  norm  that  will  be  defined  below. 

Definition  2.1: 

Let  Z  be  a  Banach  space.  Then  we  denote  by  Hn(Z)  the  space  of  all 
functions  R  +  Z  whose  first  n  derivatives  are  square-integrable  in  the 
Bochner  sense.  Moreover,  let 

Y°'n(Z)  =  {u  s  R  Z  |  e0tu,  e"0tu  e  Hn(R,  Z)  > 

and 

Xff'n(Z)  *  (u  s  R  ♦  Z  |  e  0tu  e  Hn(R,Z),  au,  e  Z  such  that 
eat(u-uso)  e  Hn(a,z)} 

The  natural  norms  in  Y°'n(Z)  and  Xff,n(Z)  are 

lul  =  le^u*  +  le  0tul 

Y  ,n(Z)  Hn(*,Z)  Hn(R,Z) 

and 

»u«  =  le  0tu»  +  lu  I  +  le0t(u-u  )• 

xa,n(Z)  Hn(R,Z)  Hn(R,Z) 

The  boundary  problem  has  been  discussed  both  analytically  and  numerically 

in  an  earlier  paper  [5] .  It  was  shown  there  that,  if  or  is  chosen  small 

enough,  then,  for  any  given  a  e  R+  and  f  e  Y°'n(H)  with  sufficiently  small 

o  #n 

norm,  there  is  a  unique  solution  u^  satisfying  u^  -  a  e  X  (R) .  Moreover 

•u^  -  al  converges  to  zero  as  *f*  converges  to  zero.  (In  [5],  we  used 
°°  2 

L  -based  rather  than  L  -based  spaces,  but  this  does  not  affect  the 
proofs)  .  Moreover,  a  numerical  scheme  was  analyzed  and  convergence  was 
proved.  In  view  of  these  results,  we  shall  henceforth  assume  that  (1.2)  has 
been  solved,  and  examine  (1.1)  with  given  Neumann  boundary  conditions 


(2.1)  u  *  b„(t)  at  x  =  +1 ,  u  *  b  „{t)  at  x  ■  -1 

x  1  x  -1 


where  b^  “  0l1 »  b_1  -  a_^  e  xa'n(R)  for  certain  ,  a  e  R+. 

This  problem  was  analyzed  in  chapters  3  and  4  of  [7] ,  and  in  the 
following  we  quote  those  results  needed  in  the  present  paper.  Following  [7], 
we  transform  (1.1)  to  a  system  of  equations,  which  can  be  classified  as 
quasi linear  parabolic.  For  this  we  substitute 


(2.2) 


p  =  u 
q  =*  u 
r  =  u 


x 

XX 


-X (t-s) 
e 

-X^t-s) 

e 

-X^t-s) 

e 

-X^t-s) 
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(u  (a)  -  u  (t) )ds 

X  X 


u  (s) 

f  XX 

ux3(s) 


Uxx(t)ux(sN 

ux4(t) 


(u"2(s) 


ux2(t))ds 


(u  (s) 

v  XX 


u  (t)u  (t) 

XX  X 

»>> 


ds 


•)ds 


With  these  substitutions,  (1.1)  is  transformed  into  the  following  system 
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p  “  r„ 


q  -  r. 


•1-2  -3  ri1?i 

r  -  -(3Hp  r  -  6Hp  qr  -  2p  l  Kg  -  -r  £  *,94) 

P  xx  x  i-l  1  ^  i-l  1  4 


(2.3) 


•i  .  i  rx 
9i  ■  -xi9i  - 


P  ip  i 


•i  \  i  4.  2rx 

g3  "  _Xig3  +  T  3 


•i  ,  i  2(  i  A  1 _ ■>  0  (  i  ^  1 _ 1 

g4  "  -Xig4  '  rxxP  93  +  7^J  *  *V«»3  +  TTJ 


The  boundary  conditions  are  now 


(2.4) 


p  -  b^  (t) ,  rx  -  (t)  at  x  -  1 

p  »  b-1  (t) ,  rx  -  b_,j(t)  at  x  “ 


(let  us  assume  that  b^  -  0^  ,  b^  -  o_1  are  at  least  In  X  '  (R) ) .  Since 
the  boundary  condition  for  p  follows  from  (2.3)  and  the  boundary  condition 
for  rx  ,  provided  it  is  satisfied  initially,  we  can  ignore  it. 

If  b  *  0,  b  .  *  0,  then,  for  any  p  e  H2[-1,1]  with  p  >  0,  a 
1-1  0  3p  0 

trivial  solution  of  (2.3)  is  given  by  p  *  pQ,  q  ■  r  «  0, 

g*  *  g*  ■  g*  *  g*  ■  0.  For  the  following,  it  is  convenient  to  set 

•  •  •  • 

*  bi  “  b-i  2  b1  +  b-1 

r  -  - - — — -  x  +  - - -  x,  so  that  the  boundary  condition  for 

4  2 

A  A 

r  is  r  *  0  at  x  ■  ±1 .  Moreover,  let  y  denote 
x 


H 


C4T 

\ 


ii 


(P_P0*  q~q0  »r •  By  A  we  denote  the  linearization  of  the  right 
hand  side  of  (2.3)  at  the  trivial  solution,  i.e.  the  operator 


i  i  i  i,  .  3h  -2a  6n  -3  A  2p0  r 


X  .1  KiV-vi-^  -v 


v  i=i 


r 

XX 

p  P0  ^0  x  p 

L  i 
i=1 

A 

A 

A 

r 

XX 

4rxg0  .  i 

2r 

x 

n  A 

p0  i 

*  4.  '  "VS  * 

poxi 

,  3 
A.  p. 
l  0 

,  i  rxx  2rxq0-| 

-  Xig4  "  T  "  pTI  ' 
1  *0  1 


The  boundary  r^  =  0  is  incorporated  in  A.  Let  denote  the  space 

H2[-1,1]  x  (H  ^  [-1 , 1  ] ) 2  x  (H1 [-1 , 1 ] )4M,  and  let  D(A),  N (A)  and  R(A)  denote 
he  domain,  nullspace  and  range  of  A,  regarded  as  an  operator  in  .  We 
quote  the  following  result  from  [7]  . 

Theorem  2.1; 

2  3p0 
Let  Pg  e  H  [-1,1]  be  given  such  that  pQ  >  0,  and  let  qQ  =  -g^-. 

Moreover,  let  n  be  an  integer  >  1.  Then,  if  o  is  chosen  small  enough  and 

•  •  o,n 

b^ (t)  and  b  (t)  have  sufficiently  small  norm  in  Y  (R) ,  equations 

(2.3),  (2.4)  have  a  solution,  for  which  y  (as  defined  above)  is 

in  X  ,n  (N  ( A) )  ©  (Y  ,n(D(A)  OR(A))  ^  Y*7'0  ^  (R(  A)  )  •  y  depends  smoothly  on 

b. 

Remarks : 

1.  If  b  is  sufficiently  smooth,  it  is  easy  to  use  a  bootstrapping  argument 
to  estimate  higher  derivatives  of  y  with  respect  to  x.  We  shall  not 
demonstrate  this  here,  but  we  shall  use  it  in  the  consistency  proof  for 
the  numerical  method. 

2.  In  [7],  the  above  theorem  was  only  formulated  for  the  special  case 

Pg  =  1  .  The  proof,  however,  carries  over  to  the  general  case. 
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The  "rubber like  liquid"  constitutive  relation  has  not  been  entirely 
sucessful  in  comparison  with  experiments,  and  several  modifications  have 
been  suggested  (see  e.g.  [3]).  Some  of  these  modified  equations  have  a 
similar  mathematical  structure  and  are  accessible  to  the  same  kind  of 
analysis.  In  this  paper  we  confine  our  attention  to  the  rubberlike  liquid 
as  the  simplest  case . 


3 .  THE  DISCRETE  PROBLEM 


In  this  chapter  we  set  up  a  discretization  scheme  for  (1.1),  (1.2)  and 

prove  convergence  for  small  forces  f  and  the  initial  condition  u(x,-°°)  =  x. 

The  infinite  interval,  on  which  (1.1)  -  (1.3)  is  posed,  is  out  at  a  finite 

point  -T  <<  0,  and  the  solution  on  (-“,T]  is  approximated  by  its  initial 

T 

value.  We  thus  get  an  approximation  u  as  solution  of  the  following 
equation : 


(3.1) 


Pu 


_  -2  -T  a  uT  (t) 

"T  »  ik(-  V  *  •'fe'TTTT^ 

u  (0  (x)) 


~  V-  /  a  ( t-s )  ( 


-T 


T 

Ux(t) 

(U^(s))2 


T 

UX(3) 

(U^(t))2 


)ds  =  0 


Jg-M.  ) 

T  2J 

(u  (t)r 

X 


for  ( x , t )  e  [-1,1]  x  [-T ,") 


(3.2)  uT(x,-T)  =  ^(x) 


-T  (uT(t))3 

(3.3)  3hu  +  /  a(t-s)ds  *  ( — 7 - -  -  V  (x)) 

X  (*  (x))Z 


r  i-  x  T  T  2 

+  /  a  ( t-s  )  [ — - -  -u(s))-f(t)(u(t))  =0 

-T  (u  (s))  X 

x 


at  x  =  +1 ,  t  >  -T. 

As  a  next  step,  this  "finite  interval  problem"  will  be  approximated  by  a 
finite  difference  procedure.  In  [5],  we  studied  the  boundary  problem  (3.3), 
using  an  implicit  Euler  scheme.  Convergence  uniformly  in  n  was  proved, 
moreover,  it  was  shown  that  the  discretization  preserves  the  property  that, 
for  appropriate  f’s,  solutions  converge  to  limits  exponentially  as  t  +  ±°®. 
Here,  we  shall  use  the  same  kind  of  scheme  for  the  full  spatio-temporal 
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problem  (3.1)  -  (3.3).  Let  us  first  introduce  some  notation 


Definition  3.1: 

Let  N  e  N,  h  ”  n'  and  for  n  e  z  let  u°  =  . VVh* 

e  R2N+3.  Then  we  define  the  "spatial"  difference  quotients 


n  n 


h 

n 

n 

ui  “ 

ui-i 

h 

n 

n 

ui+1 

-  ui-i 

2h 


and  the  "temporal"  difference  quotients 


n+1  n 
.+  n  Ui  “  ui 
6ui  — k — 


n  n-1 
n  Ui  '  Ui 
5ui  “ - k - 


n+1  n-1 

n  ui  ~  ui 
*i  *  2k 


where  k  >  0  is  given . 

Me  set  t  =nk,  n  e  Z,  t  “  -T 
n  -m 

(3.1)  is  then  discetized  as  follows  (u^ 

u(x  ,t  ) ) . 
i  n 


and  =  ih,  i  *  (-N-1 ) ( 1 ) (N+1 ) . 
denotes  the  approximation  to 
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..**«*«* 


p6+5  u?  -  3n<5+A+(-  - 

1  A"..n 


'(-4-j)  -  k  I 


(3.4) 


.  A"u"+1  • (x. )  n+1 

.  A  (-,  - - .  n|y2)  -  k  l  a(tn+rtj) 

v>  (xi)2  (A  uj  V  j— 0+1  n  3 

A~  n+l  A~  j 

A  u ,  A  u , 


a  u  A  u' 

A  CT?17  ~  ,A~  n+1 .2^  "  ° 
(A  u')  (A  u 


for  i=-N(1)N  n  >  -m. 


We  need  two  initial  vectors 


(3.5)  =  u[m_1  =  <fi(x±),  1  -  (-N-1)(1)(N+1) 


In  3.4  we  approximated  /  a(t-s)ds  by  the  discrete  sum,  although  we 
could  evaluate  it  exactly.  The  reason  is  that  this  facilitates  the  stability 


analysis  of  the  schem.  The  sum  is  given  by 


(3.6) 


k  l  a<U,  '  t .)  *  I 


M  Kflke 

r  X 


-Vtn+rt-„> 


n+1  j 


1  -  e 


Let  y",  y°1  denote  an  approximation  to  the  boundary  oondition 
uj(1).  uVl).  We  discretize  the  boundary  problem  (3.3)  as  in  [5] 


n+1  n 

y»  -  y » 


(3.7) 


k  I  -  vC—i — t  -  *  <*>) 


/  n+M 

(yA  ) 


k  j  =  -»  n+1  3  (,  (*)2 

,  n+1 .3 

n+1  (y «  )  a.  , 

+  k  l  «(t  -t  )(-— 3-7-  -  yt)  -  «tn+1)<y£  >  -  o 

j=m+1  J  (y£) 


for  £  =  +1,  -1,  n  *  -m  • 
The  corresponding  initial  conditions  are 
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(3.8) 


.  „  .  —Ml 

y1  ■  'f  (i) ,  y_1  =*  ‘fi  (-D 


Since  we  want  a  second  order  method  with  respect  to  k,  we  approximate  the 
Neumann  conditions  by  the  symmetric  difference  quotient 


(3.9) 


n+1  n+1  .  n+1  n+1 

'  y-1  '  AUN  "  yi 


The  exterior  grid  points  (-1-h,  t  1  and  (1+h,  t  „ )  can  be  eliminated  by 

n+1  n+1 

solving  (3.9)  for  u”*^  and  u^*|  and  inserting  these  values  into  (3.4). 

An  efficient  way  to  evaluate  the  last  sum  in  (3.4)  (the  discrete 
approximation  to  the  integral)  is  as  follows.  Let 


(3.10) 


(3.11) 


n+1  „  .  f 

Hi  =  K^k  l 

j— ra+1 


-VvrV  i 


(A  u*)2 


.n+1  v  v  T  A)^(tn+1  tj>A-  j 
h£  -  Ktk  l  e  A  uj 

j=-m+1 


Then  the  following  recursions  hold 


(3*12)  H+i  •  Hi  -  V  e  *  +  (e  1  *  1)gu  *  gu+1  =  0 

(A  u ) 


(3. ,3)  -  hjt  -  v  .  ‘‘V.J  ♦  (•  X,k  -  i)h;x  . 


and  since  the  n+1-st  summand  cancels  out,  (3.4)  takes  the  form 


,  -  -m 

(3.14)  p  5Tu”  -  30  5 Y(-  -L.)  -k  l  a(t  .  -t  ) 
1  A"un  j—  n+1  J 


A'u”+1  ^  (x  )  M  hj+1 

A+  (_r— i - i - )  -  A  +  (  [  g"+1A“un+1 _ ii _ ' 

Voc/  (A-u^1,2  U 1  ?Ai  1  ( A”u^+1  )2 
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a 


The  boundary  problem  can  be  treated  in  an  analogous  way.  Thus  the 
computational  solution  proceeds  as  follows.  At  a  particular  time  step,  say 
n+1,  yn+1  is  calculated  from  (3.7)  (see  [5]).  Then  g^*1  ,  h^*1  are 
calculated  from  (3.12),  (3.13)  using  the  previous  values  (time  step  n)  .  By 
inserting  this  into  (3.14),  a  nonlinear  (2N+1 )-dimensional  system  of 
equations  is  obtained  (after  substituting  u^+^ ,  according  to  (3.9)). 

This  nonlinear  system  is  solved  by  Newton's  method,  using  the  values  from  the 
previous  time  step  as  a  starting  point  for  the  iteration .  In  each  iteration 
step  the  Jacobian  is  tridiagonal  and  diagonal  dominant  (the  Gauss  algotithm 
does  not  require  pivoting) . 

In  the  remainder  of  this  chapter,  we  restrict  ourselveB  to  the  initial 
condition  'fi(x)  •  x.  We  presented  the  general  problem,  however,  since  we  have 
also  done  computations  with  different  initial  values.  For  v(x)  ■  x,  u  is 
odd,  and  the  procedure  can  be  simplified  by  imposing  the  Dirichlet  condition. 


(3.15) 


u  *  0 


at  x  =  0 


and  by  solving  in  [0,1]  x  R  .  This  cuts  the  dimension  of  the  system  down 
to  N+1.  We  shall  give  a  convergence  analysis  for  this  case. 

Since  we  also  want  to  approximate  derivatives  of  u,  we  work  in  discrete 
Sobolev  spaces  .  We  set 


(3.16)  \  *  <o0 . Vl>  e  R  lu0  '  °-  Vi  '  Vi1 


lul  -  (h  £  lu ,f)  '2 

LT  i*1 


(3.17) 


*  (u  e  ti^} ,  Ku*  2  =  2  +  2  +  *A+A  ul 


where  Au,  A+A  u  are  taken  componentwise  and  AutJ.  =  A+A~u  =  0  has  been 

N+i  N+1 
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set 


He  also  define  for  an  arbitrary  Banach  space  X: 


(3.18)  L^(X)  -  {u  -  (u")"^  I  u"  e  X  ,  lul  -  (k  l  luV)  1/2  <  •> 

n"  l£(X)  n»-*» 


and  Inductively 


(3.19)  h£(x)  -  (u  -  (uV_w  |un  e  x,  U  e  H^fx),  <S+u  e  hJ‘\x)}, 


lul  -  lul 

h£(X)  -n 


+  ld+ul 

Vi(x»  <-i  <X> 


where  H°(X)  -  L2(X) . 

Our  analysis  will  proceed  in  the  following  exponentially  weighted  Sobolev 
spaces : 


0*20}  ,»  -  <»  -  <u"f  ..l.e^vr..  e  n  •  < 


.2.2.  n  „1  ...2.  _  •  ..2 


Cnk  n  •  2  2^12, 

sucii  that  (e  (u  -  u  ))  _  e  (C(L  )  n  H.  (H.  )} 

n=»  -  k  n  k  n 


lul 


-onk  n”  ,  ,  . ,  -onk  n“  , 

x“  n=  h2(l2)  n'  hSh2* 

Xk,h  VV  VV 


+  l»"l  2  +  Me°nk(un  -  u”)f_J  +  Keank(un  -  u*))^l 

S  "X>  »v(»h  > 


(3.21) 


f  .  n. "  -onk  n  *  ,  ank  n,°*  _  ,2.2.  , 

Vk,h  *  (“  *  <u  ’n—  l<e  “  ’n - <*  “  e  W1’ 


■  ■  a,  ank  nk®  1  ^  it  ~°n*  n*  , 

M  -  ■(.  u  I  k  K.  »  2.2. 


k,h 


w 


w 
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(3.22)z£  =  { {f  =  (fn)”=_J(e"0nkfn)  e  nj(«),3f"e  R3(e0nk(fn  -  f"))  e  H^(R)} 


«f«  _  =  »(e'°nkfn)l  +  l(eank(fn  -  f"))l  +  |f  | 

Zk  «*<■» 


These  discrete  spaces  are  defined  In  an  analogous  way  as  the  spaces  used  for 
the  analysis  of  the  continuous  problem  [7]  . 

We  set  v"  =  u”  -  x^yn,  which  transforms  (3.4)  to  the  equation 

A 

F  ,  (v)  =  0,  where 

h,k,y 


(3.23)  (a)  (F  (v))”+1  =  p5+6  v"  +  px  i+i  y°  -  3n$+A+( - - - ) 

h»k»y  i  i  i  A  v  ?y  " 


.  n+1  n+1 

n+1  ,  A  v.  +  y 

-  k  l  A  : 


A  V{  +  y^ 


j- 


n+1  ~j '  '‘/a-  3  .  jv2  ,a“  n+1  ■  n+1. 2 

(A  v;:  +  y  )  (A 


\  *  y  ) 


for  i  -  1(1 )N 


(b)  (Fh,k,y(v)C  -  Vk#y<;C! 


(c)  (F,  .  (v))"+1  «  0 

h,k,y  0 


This  formula  holds  also  for  n+1  <  -m,  where  we  have  yn+1  “  1 ,  v"+1  *  0. 

A  simple  but  rather  lengthy  calculation  shows  that 


(3.24) 


h,k,y 


o’  o’ 

K  v.  *  u 

k,h  k,h 


n  *  o' 

provided  that  (y  )  _„  6  Z^  .  Considerations  very  similar  to  those  in  [5] 

<3 . 1 

show  that  this  is  the  case  if  f  (as  of  (3.3))  is  in  Y  (R) ,  a'  <  a  and 
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1 


•f •  is  small  enough.  The  boundary  conditions  at  x  *  0,1  are  already 

Y 

built  into  the  spaces.  The  main  tool  in  our  convergence  proof  is  Keller's 

nonlinear  stability-consistency  concept. 

Consistency  follows  from  a  rather  lengthy  but  essential  trivial 

calculation,  which  we  omit.  Let  v(x. ,t  )  denote  u(x, ,t  )  -  x.u  (1,t  ), 

in  in  i  x  n 

■v 

where  u  is  the  solution  of  (1.1)-(1.3)  with  <fi  (x)  *  x,  and  let  v  be  the 

gridfunction  with  values  v(x. ,t  )  in  (x. ,t  ).  Then  the  following  estimate 

in  in 


can  be  obtained. 


-(o-c') |t 


(3.25)  Ip.  .  (v)  I  <  const.  Oi  +  k  +  o(e  ~  )) 

h,k,y  o' 

Yk,h 

provided  u  -  xux(1,t)  e  xa,3(L2)  n  Xff,1(H4).  This  holds  for  f  e  Y°,4(*) 
with  sufficiently  small  norm.  The  constant  in  (3.25)  is  independent  of  h,  k, 
t-n,/  P»  and  o'  e  [Cg,^]  C  (0,<J).  Independence  of  n  is  not  assured,  since 
we  do  not  know  anything  about  the  behaviour  of  the  solution  as  n  "►  0. 

In  deriving  (3.25),  the  error  estimates  obtained  in  [5]  for  the 
approximation  (3.7),  (3.8)  of  the  boundary  condition  has  been  used  (the  choice 

CO 

of  Sobolev  rather  than  L  -based  norms  does  not  affect  the  validity  of  the 
reasoning  given  there) .  As  far  as  the  local  discretization  error  is 
concerned,  our  scheme  is  thus  second  order  accurate  in  space  and  first  order 
accurate  in  time . 

A 

In  order  to  investigate  stability,  we  linearize  F  at  v  =  0 

h,k,y  0 

(corresponding  to  u  =  x  or  f  *  0  respectively) .  We  set 


(3.26) 


and  investigate  the  equation 


^h.k.l  Fh,k, 


SV  .(0) 


(3.27) 


this  yields 


\,k,1W  =  g'  *  =  (gi>  e  Yk,h'  w  =  (wi) 
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r 


7 


(3.28)  (a) 


p5  5  w"  *  3n«  A  A  w”  +  3k  £  a(t  -  t . ) 
1  1  j— •  n+1  J 


.*+*■  n+1  a+a"  j,  ^  nH 

(A  A  wA  -  A  A  wp  +  gi 


,,  ,  n+1  .  a 

(b)  wQ  *  0,  Aw^  =  0 


we  transform  (3.28)  to  a  "parabolic"  system  of  difference  equations  by 


substituting 


,  n  -A.  (t  -  t . )  . 

n  a+a“  n  n  n  n  v  *  n  j  j 

(3.29)  q  -  A  A  w^,  r  -  6  w^  pAi=  kK^  I  e  q^ 

j=-“ 


for  i  -  1 ( 1 )M 


Then  (3.28)  assumes  the  form 


t 


(3.30)  (a)  P«+r"  -  30AVr?+1  +  3k  [  a(t 

1  1  j=-»  n+1 

-A^k 

...  jc+  n  „  n+1  .  e  -In 

(b)  6  pu  =  K*  q±  + - - -  pu 


,  .  .+  n  a+a-  n+1 

(c)  o  =*  A  A  rA 


,  n+i  -  V  n+1  n+1 

i,qi  - 3  J,  p«i  +  9i 


...  *  n+1  „  n+1 

(d)  Ar  =  0,  r.  =  0 
N  U 


The  system  (3.30)  is  reduced  further  by  a  discrete  "separation  of  variables", 
using  the  Fourier  method  (see  [8]). 


We  set 


n-  1 

(3.31)  wn  =  l  w"sin  ((2V+1)  -  x  ) 

1  v=0  Z  1 

and  likewise  for  r^,  q^,  g^,  (with  obvious  notations  etc.).  The 


•  ►-  » .*i 


-17- 


boundary  conditions  are  automatically  satisfied  by  these  expansions 
With 


(3.32)  -  (r",  q”,  R1^,...,  pjj  ),  G"  -  (g",  0,  0,...,  0)Tf 


the  system  (3.30)  can  be  rewritten  as 


(3.33) 


Av(h,k)Z^+1  -  z"  +  kG^+\  V«0(1)(N-1), 


where  Ay  is  the  (M+2) -square  matrix 


(3.34)  Av(h,k) 


1  _  /v\ 

-  kay(h) 


,.2  n+1 

—  I 

P  L 


a(t 


■}=-« 


-kX  e 


n+1 


V 


■kv 


XMk 


t  )  — 

j  P 


V 


3k 

P 


XMk 


Here  we  have  put 


(3.35) 


Vh)  = 


2V+1 

2[cos((~-i)ir)  -  1] 


n+1 

The  reader  is  reminded  that  \  *^^+1  ~  fc j )  and  hence  Ay(h,k)  does  not 

depend  on  n . 

The  stability  proof  relies  crucially  on  the  following  result. 
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Lemma  3.1: 


Av(h,k)  -  I 

The  matrix  - — -  has  M  +  1  eigenvalues  X^(h,k,v)  with  positive 

real  part  (with  a  lower  bound  for  the  real  part  not  depending  on  h,  k,  or  v) 

and  a  simple  eigenvalue  zero.  Moreover,  there  exists  a  matrix  Ev(h,k),  which 

together  with  its  inverse  is  uniformly  bounded  for  h,k  &  (0,k  ),  p  e  (0,P  ], 

Ay(h,k)  -I 

v  *=  0(1)  (N-1),  such  that  Ey(h,k).  - ^ -  Ey  (h,k)  has  the  form 


Y 

0 

0 

0 


c(*.  P.  3^7) 


where  Y  is  positive  with  a  lower  bound  uniform  in  h,?:,P  and  v  and  C  is 
analytic  in  all  its  arguments. 

Proof : 

Ay(h,k)  -  I 

The  characteristic  equation  of  - r -  is  given  by 


(3.36) 


M  K  e 

,<*>  =-f  vhl  ( ! .  -nr- 


V 


i-1  ir 

e  -  1 


n+1 

-  k  E  •(tn+l"tj,) 

jas— 00  J 


-  z 


+  Vh)  +  z)  *  0 


we  note  that  u^(h)  <  -c  for  some  c  >  0,  which  does  not  depend  on  v  or 

h.  For  convenience,  we  assume  that  X  <  X  <  ...  <  X  .  obviously 

1  2  n 

xik  _  _ 

p  (£ - - —  +  )  =  +»,  p  (2 — -  -)  =  Therefore,  there  is  at  least  one 

v  k  ^  k 
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zero  of  py  in  ( 


V 

e  -  1 


V,k  , 

e  -  1 


)  for  i  *  1(1)(M-1).  Moreover, 


pv(0)  ■  0,  and  z  =  0  is  a  simple  zero.  We  have  thus  accounted  for  M 
zeros  of  py,  and  two  more  are  left. 

We  have 


(3.37) 


Pv(z)’fVh)  % 


M 

I 

i-1 


V 


V 


V 

(*— -  z)2 


+  f-  «w(h)  +  2z 


For  z  e  R  and  z  <  0,  we  find  py(z)  <  0,  whence  py  has  no  negative 
roots .  On  the  other  hand,  we  find 


(3 .38) Ira  pv(z) 


-  In  z  (ay(h) 


M 

I 

i=1 


V 

V 

- 

,e  1  -  1 


+  n  av(h)  +  Re  z) 


-  z 


If  Im  py(z)  =*  0,  we  therefore  have  either  Im  z  *  0  or  Re  z  >  c  >  0.  if  z 
is  real  and  py(z)  *  0,  then  we  know  z  cannot  be  negative,  and  if  z  is 
positive,  one  easily  finds  a  lower  bound  for  Re  z  from  (3.36).  Thus  we  have 
proved  the  first  part  of  the  lemma. 

The  second  part  follows  from  elementary  matrix  manipulations,  which  we 

3n  •  • 

don  not  demonstrate  here.  One  finds  that  Y  =  -  a  (h)  +  Y  where  Y  is 

P  v 

analytic  in  k,  p  and  — rr-r  • 

Vh) 

Lemma  3.2* 

Let  the  nxn-matrix  C  have  only  eigenvalues  with  positive  real  part. 
Then,  for  each  a  €  L^< (-“,")?  Rn),  the  difference  equation 

(3.39)  <S"yit)  -  y(--)  =•  -Cy( t)  +  a(t) 


20 


2  n 

has  a  unique  solution  y  e  L  ((-*,0®);  R  )  and  the  estimate 

(3.40)  lyl  <  c«a« 

L 2  L2 

holds  with  c  >  0  independent  of  k  and  a. 

Proofs 

Let  y  and  a  be  the  Fourier  transforms  of  y  and  a: 


1  f  x  2irist.  ....  1  f  x  2irist 

y(t)  ■ -  J  y(s)e  ds,  a(t)  »  -  J  a(s)e  ds 

/2n  -*>  /27T 


Equation  (3.39)  is  equivalent  to 

-2irik 

1  -  e 


y(s)  »  -Cy(s)  +  a(s) 


or 


-2*ik  _  _ 

y(s)  =  (  - -  I  +  c)  a(s) 


Since  Re ( 


1  -  e 


-2irik 


-)  >  0,  the  lemma  is  immediate  (use  Parseval's  equality). 


As  a  consequence,  we  have  15  yl  <  const.  Ia<  . 

L2  L2 


An  easy  corollary  is  that  for  a  e  Y°'^(Rn),  there  is  a  unique  solution 
y  e  Y°'°(Rn)  and  an  estimate  *y ■  <  dal  holds,  provided  a  is 

*  *  tal' 

small  enough  (to  see  this,  substitute  y  *  e  y  and  apply  lemma  3.2).  If 

0,0  n  o,0  n 

C  -  0  and  a  e  Y  (R  ) ,  there  is  a  solution  y  e  X  (R  )  and  an  estimate 

lyl  <  dal  holds  (recall  definition  2.1  for  Y0,^,X°'^).  From  these 

X  '  Y  * 

results  the  estimate 


(3.41) 


r+  n. 


I (Z  )  I  .  +  15  z  I 


< 

is  immediate.  The  constant  is  independent  of  h,  k,  v,  and  P,  but  may 


-2  1  - 


depend  on  n.  When  we  rewrite  the  norms  in  x°  . 

k,h 


in  terms  of 


o  ’ 

and  Y,  . 
k,h 

Fourier  coefficients,  we  immediately  see  that  the  solution  of  (3.27)  satisfies 
an  estimate 

A  A 

(3.42)  Iwl  0,  <  const.  Igi  ,  , 

k,h  k,h 


i.e.  we  have  stability  unconditionally  with  respect  to  k  and  h,  and  the 
constant  is  uniform  in  h,  k,  p. 

A  simple  calculation  shows  that,  in  a  neighbourhood  of  v  *  0,  y  *  1,  the 

•  A 

derivative  F^  (v)  is  uniformly  Lipschitz  continuous  with  respect  to  y 

A 

and  v,  i.e.  we  have 


(3.43) 


y 


(v) 


Fh.k.y(v*V 

k,h 


♦  Y 


o' 

k,h 


<  const.  ( lv 


O' 

k,h 


■y  -  y*'  0.) 

Z, 


If 

(3.7), 


v 

(3.8) 


denotes  the  solution  to  the  exact  problem,  and 
>■  ^ 

,  then  ly  -  ll  _  and  Ivl  .  are  small  if 
0  o' 

*k  Xk,h 


y 

If! 


the  solution  to 

0,2  18 


small.  Therefore  all  conditions  for  Keller's  [1]  theory  are  satisfied,  and  we 
obtain  the  final  result. 

Theorem  3.1: 

Let  o  be  a  positive  real  number  chosen  sufficiently  small.  For  any 
f  e  Y°'4(R),  which  has  sufficiently  small  norm  in  Y°'2(R),  there  is  a  unique 
solution  ( (u8) )  _* )  of  (3.4),  (3.5),  (3.9),  and  the  convergence 

estimate 


(3.44) 


■(«•:- «<v 

Xk,h 

2  -(0-0*) |t_m| 

<  const,  (h  ♦  k  +  o(e  ) 


-22- 


\ 


holds  as  h  +  0,  k  ♦  0,  t  The  constant  is  independent  of 

-in 

h,  k,  o'  e  [0^,0^]  c  (0,0).  Moreover,  (3.4),  (3.5),  (3.9)  can  be  solved  by 

the  Newton  procedure  (going  from  the  nth  to  the  (n+1)st  time  step),  which 

is  quadratically  convergent  from  a  sphere  of  starting  values  which  does  not 

shrink  to  a  point  as  h  -►  0,  k  +  0,  t  ♦  t  *  +“.  The  convergence  is 

-tn  n 

unconditional  in  the  sense  that  h,  k,  t_m  may  approach  0  (or  -*• 
respectively)  independently. 

The  main  restriction  of  this  theorem  is  the  required  smallness  of  f . 
However  not  more  can  be  expected  theoretically  since  the  existence  theory  only 
holds  for  small  f's.  In  the  next  chapter  we  report  computation  using 
arbitrary  (actually  large)  forces  f  and  convergence  was  realized  numerically 


4.  DISCUSSION  OF  NUMERICAL  RESULTS 

For  the  computations  we  used  the  same  kernel  as  '.n  [5]  : 

8  -At 
a(t)  *  l  K  e  1 
i=1 

with  the  following  constants  ,  A^  . 


i 

Ai(sec“1 ) 

Ki(Nm“2sec“1 

1 

CO 

1 

o 

1  x  10-3 

2 

10"2 

CN 

1 

O 

X 

CO 

• 

3 

10"1 

1.89  x  10_1 

4 

1 

9.8  x  103 

5 

10 

2.67  x  105 

6 

102 

5.86  x  106 

7 

103 

9.48  x  107 

8 

104 

1  .29  x  109 

These  numbers  were  obtained  by  Laun  [2]  from  an  experimental  fit  for  a 

polyethylene  melt  at  150°  C,  which  he  calls  "Melt  1". 

The  parameter  n  is  physically  identified  as  a  Newtonian  contribution  to 

the  viscosity.  Experimental  values  are  not  available,  and  theoretically  n 

is  either  a  solvent  viscosity  (for  polymer  solutions)  or  it  results  from 

fractions  of  low  molecular  weight  (for  melts)  .  The  value  of  n  should  be 

compared  to  the  viscosity  resulting  from  the  memory,  which,  for  constant  shear 

§  _2  _2 

rate,  is  given  by  l  K  A,  38  50000Nm  sec.  One  would  expect  n  to 
i=1  1 

influence  the  solution  significantly  only  if  it  is  at  least  comparable  to  this 
value.  This  heuristic  argument  was  confirmed  in  [5]  for  the  boundary 
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equation,  and  also  by  our  present  computations  for  the  spatial  problem.  The 

numbers  given  in  the  following  are  understood  to  be  in  the  following  units: 

-2 

0  is  given  in  Nm  sec. 

f  denotes  the  force  acting  on  the  ends  of  the  filament  divided  by  the 
cross-sectional  area  in  the  reference  state  (u  *  x)  .  It  is  measured  in 
Nm  2  . 

P  denotes  the  density  of  the  filament  multiplied  by  the  square  of  half 
the  length  in  the  reference  state  (this  latter  scaling  factor  arises  from  the 
normalization  of  the  variable  x  to  the  interval  [-1,1]).  P  is  measured  in 
kg  m 

Time  is  measured  in  seconds .  The  mesh  size  for  the  following  plots  was 

-t2/25 

0.1  for  both  t  and  x.  For  plots  1-13  we  chose  f  =  100000  e  and 

u(x,-»)  =  x. 

The  first  three  plots  show  u,  u  and  u  for  p  =  1 ,  0=1.  It  can 

X  XX 

be  seen  that  Uxx  is  negligibly  small,  and  u  is  almost  linear  in  x.  That 
means  the  solution  is  determined  by  the  evolution  of  the  boundary  condition, 
and  inertial  forces  can  be  neglected. 

This  changes,  if  P  is  increased.  Physically,  this  means  changing  the 
length  of  the  filament.  For  realistic  values  of  the  density,  P  =  1  would 
correspond  to  an  initial  length  of  a  few  millimeters,  p  =  1000  would 
correspond  to  an  initial  length  of  about  1m. 

In  plots  4-6,  we  have  p  =  1000,  n=1.  Three-dimensional  plots  of  u 

and  u  and  sections  of  u  at  various  times  are  shown.  It  can  be  seen 

XX  X 

that  u^x  has  increased  by  a  factor  of  1000  compared  to  the  previous 
plots.  Otherwise  the  qualitative  behaviour  remains  roughly  the  same. 

The  next  three  plots  were  made  for  the  same  P  and  h  =  100000  (for 
0  <  0  <  10000,  the  solutions  changed  very  little).  In  comparison  to  0=1, 
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it  was  found  here  that  the  boundary  value  for  u^  increases  more  slowly  up 
to  t  =  -2  and  then  increases  rather  suddenly  around  t  «  0 .  In  this  region 
inertial  forces  become  very  important,  a  fact  manifested  in  the  plots  by  a 
rather  pronounced  spike  in  u  . 

In  figure  10  we  have  P  =  1000,  0  =  1000000.  In  this  case  the  behaviour 

becomes  almost  Newtonian,  and  there  is  hardly  any  elastic  recovery  (the 
maximal  value  for  u(x=1)  is  1.404,  the  value  at  t  =  60  is  still  1,398). 
The  dependence  of  u  on  x  is  again  almost  linear. 

Several  calculations  were  done  for  P  =  10000.  Figures  11-13  show 
u,  u  and  u  for  h  =  100000.  It  can  be  seen  that  u  becomes  rather 

X  XX  XX 

large.  When  one  looks  carefully  at  the  plots  for  u,  one  also  finds  that  a 
little  "overshoot"  occurs  in  the  relaxation:  Whereas  in  the  previous  plots 
u  decreased  monotonically  after  reaching  the  maximal  value,  this  is  no  longer 
true  here: 


t  =  6  t  =  8  t  —  10  t  =*  12 

u(x=-1)  -83.6  -22.7  -3.8  -4.3 

If  smaller  values  of  h  are  chosen,  this  "overshoot"  becomes  even  more 

pronounced.  The  mesh  size  becomes  very  crucial  here?  if  it  is  chosen  too 

coarse,  then  the  numerical  approximation  to  ux  changes  sign  in  finite 

time.  This  effect  could  not  be  reproduced  with  a  finer  mesh,  although  the 

values  agree  very  well  before  one  gets  near  this  "critical"  point. 

For  the  next  three  plots  we  have  chosen  an  "oscillating"  force:  f  *  0 

2 

for  t  <  -20  and  f  *  20000*  e  t  //625 ( 1 +cosj)  for  t  >  -20.  The  values  P  * 
1000,  n  =100000  were  chosen.  The  solution  seems  to  "follow"  the  oscillations 
with  a  certain  time  lag  (cf.  [5J )  . 
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In  plots  17-31  we  are  concerned  with  filaments  which  are  not  in  the 

reference  state  a  =  x  at  t  =  -*,  but  rather  in  an  already  deformed  state. 

2 

In  all  cases  we  have  chosen  f  =  100  000  e-t  ^25.  In  figures  17-19,  we  have 

3 

u(x,-°®)  =  x  +  x  ,  corresponding  to  a  filament  which  becomes  thinner  towards 

the  ends.  The  parameters  P  =  1000,  h  *  1  were  chosen.  Observe,  in 

particular,  that  the  recovery  of  u  shows  several  oscillations  here. 

The  last  twelve  plots  have  the  seme  f  and 
2 

u(x,  -°»)  =  x  -  — ,  representing  a  filament  that  is  thinner  towards  the  left. 
In  figures  20-25,  we  have  P  =  1 ,  h  =  1,  in  the  second  series  of  figures  (26- 
SI  ),  we  have  p  =  10000,  n  =  100000.  The  qualitative  behaviour  changes 
considerably,  and,  in  particular,  the  effect  of  inertia  is  very  important. 
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